data {
  // observations are split into four datasets based on whether
  // (1) message midpoint hits the boundary and
  // (2) message range hits the maximum size
  int<lower = 1> N_id;     // number of subjects
  int<lower = 1> N_period; // number of periods
  // observations in which midpoint does not hit boundary and range
  //   does not hit maximum value
  int<lower = 1> N_1;
  int<lower = 1, upper = N_id> id_1[N_1];
  int<lower = 1, upper = N_period> period_1[N_1];
  vector[N_1] Target_1;
  vector[N_1] Message_Midpoint_1;
  vector[N_1] log_Message_Range_Size_1;
  // observations in which midpoint does not hit boundary but range
  //   does hit maximum value
  int<lower = 0> N_2;
  int<lower = 1, upper = N_id> id_2[N_2];
  int<lower = 1, upper = N_period> period_2[N_2];
  vector[N_2] Target_2;
  vector[N_2] Message_Midpoint_2;
  vector[N_2] log_Max_Message_Range_Size_2;
  // observations in which midpoint hits upper boundary
  int<lower = 0> N_3;
  int<lower = 1, upper = N_id> id_3[N_3];
  int<lower = 1, upper = N_period> period_3[N_3];
  vector[N_3] Target_3;
  real U; // right censor point for Message_Midpoint
  // observations in which midpoint hits lower boundary
  int<lower = 0> N_4;
  int<lower = 1, upper = N_id> id_4[N_4];
  int<lower = 1, upper = N_period> period_4[N_4];
  vector[N_4] Target_4;
  real L; // left censor point for Message_Midpoint
}

parameters {
  real a0;                    // grand mean for intercept in midpoint model
  vector[N_id] a_e_id;        // (unscaled) subject-level random intercept
  real<lower = 0> a_s_id;     // scale for subject-level random intercept
  real b0;                    // grand mean for slope in midpoint model
  vector[N_id] b_e_id;        // (unscaled) subject-level random slope
  real<lower = 0> b_s_id;     // scale for subject-level random slope
  real<lower = 0> s_midpoint; // scale for error in midpoint model
  real g0;                    // grand mean for intercept in range model
  vector[N_id] g_e_id;        // (unscaled) subject-level random intercept
  real<lower = 0> g_s_id;     // scale for subject-level random intercept
  real h0;                    // grand mean for slope in range model
  vector[N_id] h_e_id;        // (unscaled) subject-level random slope
  real<lower = 0> h_s_id;     // scale for subject-level random slope
  real<lower = 0> s_logrange;    // scale for error in range model
}

transformed parameters {
  // containers for subject-level coefficients
  vector[N_id] a_id;
  vector[N_id] b_id;
  vector[N_id] g_id;
  vector[N_id] h_id;
  // containers for midpoint model linear predictors
  vector[N_1] m_midpoint_1;
  vector[N_2] m_midpoint_2;
  vector[N_3] m_midpoint_3;
  vector[N_4] m_midpoint_4;
  // containers for range model linear predictors
  vector[N_1] m_logrange_1;
  vector[N_2] m_logrange_2;
  vector[N_3] m_logrange_3;
  vector[N_4] m_logrange_4;
  a_id = a0 + a_s_id * (a_e_id - mean(a_e_id));
  b_id = b0 + b_s_id * (b_e_id - mean(b_e_id));
  g_id = g0 + g_s_id * (g_e_id - mean(g_e_id));
  h_id = h0 + h_s_id * (h_e_id - mean(h_e_id));
  for (n in 1:N_1) {
    m_midpoint_1[n] = a_id[id_1[n]] + b_id[id_1[n]] * Target_1[n];
    m_logrange_1[n] = g_id[id_1[n]] + h_id[id_1[n]] * Target_1[n];
  }
  for (n in 1:N_2) {
    m_midpoint_2[n] = a_id[id_2[n]] + b_id[id_2[n]] * Target_2[n];
    m_logrange_2[n] = g_id[id_2[n]] + h_id[id_2[n]] * Target_2[n];
  }
  for (n in 1:N_3) {
    m_midpoint_3[n] = a_id[id_3[n]] + b_id[id_3[n]] * Target_3[n];
    m_logrange_3[n] = g_id[id_3[n]] + h_id[id_3[n]] * Target_3[n];
  }
  for (n in 1:N_4) {
    m_midpoint_4[n] = a_id[id_4[n]] + b_id[id_4[n]] * Target_4[n];
    m_logrange_4[n] = g_id[id_4[n]] + h_id[id_4[n]] * Target_4[n];
  }
}

model {
  // prior
  a0 ~ normal(0, 10);
  a_s_id ~ cauchy(0, 2.5);
  a_e_id ~ normal(0, 1);
  b0 ~ normal(0, 10);
  b_s_id ~ cauchy(0, 2.5);
  b_e_id ~ normal(0, 1);
  s_midpoint ~ cauchy(0, 2.5);
  g0 ~ normal(0, 10);
  g_s_id ~ cauchy(0, 2.5);
  g_e_id ~ normal(0, 1);
  h0 ~ normal(0, 10);
  h_s_id ~ cauchy(0, 2.5);
  h_e_id ~ normal(0, 1);
  s_logrange ~ cauchy(0, 2.5);
  // posterior
  Message_Midpoint_1 ~ normal(m_midpoint_1, s_midpoint);
  Message_Midpoint_2 ~ normal(m_midpoint_2, s_midpoint);
  log_Message_Range_Size_1 ~ cauchy(m_logrange_1, s_logrange);
  // increment log probability for censoring
  target +=
    cauchy_lccdf(log_Max_Message_Range_Size_2 | m_logrange_2, s_logrange) +
    normal_lccdf(U | m_midpoint_3, s_midpoint) +
    normal_lcdf(L | m_midpoint_4, s_midpoint);
}
